DALS017-高维数据推断2-统计原理

发表于 | 分类于
| 字数统计: 8,443 | 阅读时长 ≈ 34

前言

这一部分是《Data Analysis for the life sciences》的第6章线性模型的第2小节，这一部分的主要内容涉及高维数据统计的一些原理，相应的R markdown文档可以参考作者的Github。

高维数据的p值

在前面我们了解了，当我们在分析高维数据时，p值就不再是一个很好的统计指标了。这是因为我们在同一时间同时检测了许多特征值(feature)。这种检测手段被称为多重比较( multiple comparison)，或多重检测(multiple testing)，或多重性问题(multiplicity)。在此情况下，p值就不再适用。另外，当我们同时检测多个假设问题时，仅仅基于一个小p值的阈值，例如0.01，这就很容易导假阳性。针对这种情况，我们需要定义一个新的术语来研究高通量数据。

为了处理多重比较的问题，我们广泛使用的方法就是定义一个程序(procedure，也可以说是一种算法，也可以翻译为校正等等，总之，表达的是一个意思)，然后用它来估计或控制(control)计算过程中的错误率(rate error)。我们这里所说的控制(control)的意思是说，我们会采用这个程序来保证错误率(error rate)低于某个提前设定的值。通过参数或截止值(cutoff)来进行设定的这个程序通常比较灵活，它会让我们能够控制特异性(specificity)和灵敏度(sensitivity)，这种程序的一个典型功能如下所示：

  • 计算每个基因的p值；
  • 计算出p值小于$\alpha$的所有显著性基因。

这里需要注意的是，当我们改变$\alpha$值时，会调整相应的特异性(specificity)和灵敏度(sensitivity)。

接着我们来定义错误率(error rate)，它会让我们对统计过程进行估计和控制。

错误率

在这一部分中，我们先了解到I类错误与II类错误，这两类错误分别代表假阳性(false positives)与假阴性(false negatives)。通常，特异性(specificity)与I类错误有关，灵敏性(sensitivity)与II类错误有关。

在一次高通量实验里，我们会犯第I类错误和第II类错误。我们参考了Benjamini-Hochberg的论文，做了以下表格，总结了这些错误，如下所示：

Called significant（真） Not called significant（假） Total
Null True $V$ $m_0-V$ $m_0$
Alternative True $S$ $m_1-S$ $m_1$
True $R$ $m-R$ $m$

为了说明这个表格中的内容，我们来打个比方，假设我们检测了10000个基因，这就意味着我们需要做检测的次数为$m=10000$。

那些零假设是真的基因数目（多数是我们不感兴趣的基因，零假设为真，也就是说这些基因在对照和实验组中都没有显著性差异）为$m_{0}$，那些零假设为假的基因数目为$m_{1}$，零假设为假，也就是说，替代假设(alternative hypothesis)为真（不一定严谨，反正就是说零假设为假）。通常来说，我们感兴趣的是尽可能地检测到那些替代假设为真的基因（真阳性），避免检测到那些零假设为真的基因（假阳性）。对于多数高通量实验来说，我们会假设$m_{0}$远大于$m_{1}$（这句话我的理解就是，在一次高通量实验中，没差异基因的数目$m_{0}$要大于有差异基因的数目$m_{1}$）。

例如我们检测了10000个基因，对其中约有100个基因感兴趣。这也就是说，$m_1 \leq 100$ 并且 $m_0 \geq 19,900$。

在这一章中，我们指的特征值(feature)就是我们的检测值。在遗传学中，这些特征值可以是基因(genes)，转录本(transcripts)，结合位点(binding sites)，CpG岛和SNPs。

在上面的那个表格中，$R$表示的是经过检测后，有显著性差异的特征值的数目总和，而$m-R$则表示不显著的基因数目。表格中剩下的部分表示的是一些实际上未知的重要的量。

  • $V$表示I类错误或假阳性。更具体地来说就是，$V$表示了那些零假设为真的特征值的数目。
  • $S$表示的是真阳性的数目。具体地来说就是，$S$表示替代假设为真的特征值的数目。

$m_{1}-S$表示了II类错误或假阴性，$m_{0}-V$表示真阴性。

这里需要牢记的是，当我们只做一次检测时，p就仅仅是当$V=1$，$m=m_{0}=1$时的概率。功效(power)就是当$S=1$，$m=m_{1}=1$时的概率。在这种非常简单的案例里，我们并不制作上面类似的表格，但是我们会说明一下，如何定义表格中的术语，从而帮助我们处理高维数据。

数据案例

现在看一个案例。在这个案例中我们会使用小鼠数据进行Monte Carlo模拟，从而模拟一种情况，在这种情况里，我们会检测10000种对小鼠体重无影响的减肥饲料(fad diets)。这就是说，在零假设下，这些饲料对小鼠体重没影响为真，也就是说$m=m_{0}=10000$，并且$m_{1}=0$，现在我们先进行一个样本数目为12的计算，并且我们认为p值小于$\alpha=0.05$显著，如下所示：

1
2
3
4
5
6
7
8
9
10
11
set.seed(1)
population = unlist( read.csv("femaleControlsPopulation.csv") )
alpha <- 0.05
N <- 12
m <- 10000
pvals <- replicate(m,{
  control = sample(population,N)
  treatment = sample(population,N)
  t.test(treatment,control)$p.value
})
sum(pvals < 0.05) ##This is R

结果如下所示：

1
2
> sum(pvals < 0.05) ##This is R
[1] 462

从结果我们可以看出，这个假阳性（462个）还是比较高的，这是要在多数分析中是要避免的。

下面我们来看一个更加复杂的代码，这段代码会进行人为设定10%的饮食有效，平均效应(average effect size)为$\Delta=3$盎司（约85克）。仔细研究这段代码可以帮助我们理解上面的那个表格，现在我们先来定义这样的数据，其中10%有效：

1
2
3
4
5
6
7
8
alpha <- 0.05
N <- 12
m <- 10000
p0 <- 0.90 ##10% of diets work, 90% don't
m0 <- m*p0
m1 <- m-m0
nullHypothesis <- c( rep(TRUE,m0), rep(FALSE,m1))
delta <- 3

现在我们做10000次模拟统计，每次都采用t检验，并记录下来：

1
2
3
4
5
6
7
8
9
set.seed(1)
calls <- sapply(1:m, function(i){
  control <- sample(population,N)
  treatment <- sample(population,N)
  if(!nullHypothesis[i]) treatment <- treatment + delta
  ifelse( t.test(treatment,control)$p.value < alpha,
          "Called Significant",
          "Not Called Significant")
})

由于我们知道哪些是数据是有差异的（毕竟是自己人为生成的，保存在了nullHypothesis中)，我们现在计算一下上面的表格，代码如下所示：

1
2
null_hypothesis <- factor( nullHypothesis, levels=c("TRUE","FALSE"))
table(null_hypothesis,calls)

结果如下所示：

1
2
3
4
5
6
> null_hypothesis <- factor( nullHypothesis, levels=c("TRUE","FALSE"))
> table(null_hypothesis,calls)
               calls
null_hypothesis Called Significant Not Called Significant
          TRUE                 421                   8579
          FALSE                520                    480

结果中的第1列就是$V$与$S$，需要注意的是，$V$与$S$是随机变量，如果我们再次运行这段代码以，这些值就会改变，如下所示：

1
2
3
4
5
6
7
8
9
10
11
B <- 10 ##number of simulations
VandS <- replicate(B,{
  calls <- sapply(1:m, function(i){
    control <- sample(population,N)
    treatment <- sample(population,N)
    if(!nullHypothesis[i]) treatment <- treatment + delta
    t.test(treatment,control)$p.val < alpha
  })
  cat("V =",sum(nullHypothesis & calls), "S =",sum(!nullHypothesis & calls),"\n")
  c(sum(nullHypothesis & calls),sum(!nullHypothesis & calls))
})

结果如下所示：

1
2
3
4
5
6
7
8
9
10
V = 410 S = 564 
V = 400 S = 552 
V = 366 S = 546 
V = 382 S = 553 
V = 372 S = 505 
V = 382 S = 530 
V = 381 S = 539 
V = 396 S = 554 
V = 380 S = 550 
V = 405 S = 569

针对这种情况，我们就定义了错误率(error rate)。例如，我们可以估计$V$大于0的概率。它可以解释为，在10000次检验中，我们出现I类错误的概率。在上述模拟数据中，$V$在每次模拟中都大于1，因此我们怀疑这个概率实际上就是1（这里的1就是“V大于0”这个事件的概率，换句话讲，就是，在这个检验中，必然会出现假阳性）。

当$m=1$时，这个概率就等于p值。当我们进行多重检验模拟时，我们称之为多重比较谬误(Family Wide Error Rate,FWER)，它与我们广泛使用的一个检验方法有关，即Bonferroni校正( Bonferroni Correction)。

Bonferroni校正

现在我们了解一下FWER是如何运用的，在实际的分析中，我们会选择一个程序(procedure)来保证FWER小于某个提前设定好的阈值，例如0.05，通常情况下，就使用$\alpha$来表示。

现在我们来描述一个这样的程序：拒绝所有p值小于0.01的假设。

为了说明这个目的，我们假设所有的检验都是独立的（在10000个饲料检验的实验里，这个假设是成立的；但是在一些遗传学实验里，这个假设有可能不成立，因为某些基因会共同发挥作用）。每次检验得到的p值我们用$p_1,\dots,p_{10000}$表示。这些独立的随机变量如下所示：

如果我们要模拟这个过程，代码如下所示：

1
2
3
B<-10000
minpval <- replicate(B, min(runif(10000,0,1))<0.01)
mean(minpval>=1)

结果如下所示：

1
2
> mean(minpval>=1)
[1] 1

因此，我们计算出来的FWER是1，这并非我们想要的结果。如果我们想要它低于$\alpha=0.05$，那么我们的统计就是失败的。

我们如何做才能使得一个错误出现的概率低于$\alpha$呢，我们可以使用上面的公式推导一下，通过选择更严格的截止值(cutoff），之前的是0.01，从而将我们的错误降低至至少5%的水平，如下所示：

解这个方程，我们就得到了 $1-(1-k)^{10000}=0.01 \implies k = 1-0.99^{1/10000} \approx 1e-6$

这就是我们给出的一个程序案例。这实际上就是Sikdak过程。尤其是，当我们定义一个说明，例如拒绝所有p值小于 0.000005的零假设。然后，我们知道了p值是一个随机变量，我们会使用统计理论来计算，如果我们遵循这个程序，平均会犯多少错误。更确切地讲就是，我们可以计算出这些错误的边界，也就是说，这些错误小于某些预定值的比例。

Sidak校正的一个问题是，这个校正假设所有的检验都是独立的。因此当这个假设时成立的，它只能控制FWER。百Bonferroini校正则更为通用，因为即使每个检测不独立，它也能控制FWER。。与Sidak校正一样，我们首先来看一下：

现在使用前面表格的那种表示方法为：

Bonferoni校正会设定$k=\alpha/m$，因此可以写为如下形式：

将FWER控制在0.05水平是一种非常保守的方法，现在我们使用前面计算的p值来计算一下：

1
2
3
4
5
6
7
set.seed(1)
pvals <- sapply(1:m, function(i){
  control <- sample(population,N)
  treatment <- sample(population,N)
  if(!nullHypothesis[i]) treatment <- treatment + delta
  t.test(treatment,control)$p.value
})

只需要关于p值是否小于0.05/10000即可，如下所示：

1
sum(pvals < 0.05/10000)

结果为：

1
2
> sum(pvals < 0.05/10000)
[1] 2

当使用了Bonferroni校正后，即使我们进行了10000次饮食实验，却只有发现2个假阳性的结果，Bonferroni是一种非常严格的校正。

错误发现率(FDR)

在许多情况下，要求FWER是0.05没多大意义，因为它太严格了。例如，我们常见的做法就是先进行初步的小型实验来确定小数候选基因。这种做法被称为之发现驱动的实验或项目。我们也许会寻找一个未知的致病基因，而不仅仅是对候选基因进行采用更多的样本进行后续研究。如果我们开发一个程序，例如一个10个基因列表，从中发现1到2个为重要的基因，这个实验就算非常成功。小样本实验，实现$FWER\leq0.05$的唯一途径就是使用一个空的基因列表。在前面我们已经看到，虽然只有1000种包含有效，但是最终我们只得到2个饮食有效这样一个结果，如果把样本数目降低到6，这个结果有可能就是0，如下所示：

1
2
3
4
5
6
7
8
set.seed(1)
pvals <- sapply(1:m, function(i){
  control <- sample(population,6)
  treatment <- sample(population,6)
  if(!nullHypothesis[i]) treatment <- treatment + delta
  t.test(treatment,control)$p.value
})
sum(pvals < 0.05/10000)

结果如下所示：

1
2
> sum(pvals < 0.05/10000)
[1] 0

由于我们要求$FWER\leq 0.05$，因此我们实际上就是0功效（灵敏度）。在许多方法中，这种情况称为特异性(specificity)过强(over-kill)。替代FWER的另外一种方法就是FDR，即错误发现率(false discover rate)。FDR背后的思想就是集中关注Q值，即 $Q \equiv V/R$，当$R=0,V=0$时，$Q=0$。其中，当$R=0$（没有显著性）就表明，$V=0$（没有假阳性）。因此$Q$是一个随机变量，它的范围是0到1，通过计算Q的平均值，我们可以定义一个比值(rate)，它所表示的是意思是，在显著的基因里，假阳性的基因占的比例。为了更好地理解这个概含，我们计算Q的程序要求调用所有的p值都小于0.05。

向量化代码

在R中可以使用sapply()系列函数来加快运行速度，前面已经看到了这个函数的使用方法，现在使用传统的代码来看一下Q值的计算方法，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
library(genefilter) ##rowttests is here
set.seed(1)
##Define groups to be used with rowttests
g <- factor( c(rep(0,N),rep(1,N)) )
B <- 1000 ##number of simulations
Qs <- replicate(B,{
  ##matrix with control data (rows are tests, columns are mice)
  controls <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
  ##matrix with control data (rows are tests, columns are mice)
  treatments <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
  ##add effect to 10% of them
  treatments[which(!nullHypothesis),]<-treatments[which(!nullHypothesis),]+delta
  ##combine to form one matrix
  dat <- cbind(controls,treatments)
  calls <- rowttests(dat,g)$p.value < alpha
  R=sum(calls)
  Q=ifelse(R>0,sum(nullHypothesis & calls)/R,0)
  return(Q)
})

结果如下所示：

1
2
3
4
> head(Qs)
[1] 0.4513274 0.4063786 0.4568627 0.4490414 0.4468314 0.4315569
> tail(Qs)
[1] 0.4390756 0.4718657 0.4395373 0.4425711 0.4536391 0.4284284

控制FDR

上述代码使用Monte Carlo模拟计算了1000次10000个样本的实验，并保存了Q值，现在看一下Q值的直方图，如下所示：

1
2
3
library(rafalib)
mypar(1,1)
hist(Qs) ##Q is a random variable, this is its distribution

FDR就是Q值的平均值，如下所示：

1
2
FDR=mean(Qs)
print(FDR)

如下所示：

1
2
> print(FDR)
[1] 0.4463354

这里的计算结果表明，FDR值比较高。这是因为对于90%的统计检验而言，零假设是真。这就暗示了，由于p值是cutoff值为0.05，100个检验中，大概有4到5个是假阳性。再加上我们没有考虑到替代所设为真时的所有情况，因此FDR值就比较高。那么我们如何控制它呢，如果我们想要更低的FDR值，比如5%怎么办？

为了用图形说明FDR为什么这么高，我们来绘制p值的直方图。我们使用一个较大的m值从直方图中获取更多的数据。再绘制一条水平线，表示当NULL为真时，从$m_{0}$个案例（指的基因或其它检测值，$m_{0}$就是没有统计学差异的基因）中到的阳性结果的均匀分布，如下所示：

1
2
3
4
5
6
7
8
9
set.seed(1)
controls <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
treatments <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
treatments[which(!nullHypothesis),]<-treatments[which(!nullHypothesis),]+delta
dat <- cbind(controls,treatments)
pvals <- rowttests(dat,g)$p.value
h <- hist(pvals,breaks=seq(0,1,0.05))
polygon(c(0,0.05,0.05,0),c(0,0,h$counts[1],h$counts[1]),col="grey")
abline(h=m0/20)

第1个柱子（灰色）表示提p值小于0.05的基因的数目，从水平线可知，大概有一半p值小于0.05的基因是假阳性，这与前面提到的FDR=0.5是一致的。我们来看一下当柱子的宽度为0.01时FDR的值会更低，但同时我们的显著性差异数目也会降低。

1
2
3
h <- hist(pvals,breaks=seq(0,1,0.01))
polygon(c(0,0.01,0.01,0),c(0,0,h$counts[1],h$counts[1]),col="grey")
abline(h=m0/100)

当柱子的宽度变为0.01时，我们可以看到一个更小的p值cutoff，并且检测到到特征值的数目也会下降（降低了灵敏性sensitivity)，但是，FDR也会降低（提高了特异性specificity）。我们如何设定这个cutoff呢？其中一个方法就是设定一个FDR水平$\alpha$，然后设置控制错误率的程序即可：$FDR \leq \alpha$。

Benjamini-Hochberg

对于任意给定的$\alpha$，Benjamini-Hochberg方法都非常适用，这种方法可以让使用者地每个检验的p值进行校正，也能很好地定义一个程序。

Benjamini-Hochberg方法的原理是，先按照升序对p值进行排列，即$p_{(1)},\dots,p_{(m)}$，然后定义一个$k$用来表示最大的秩$i$，它所对应的p值为：

这个程序会拒绝那些p值小于或等于$p_{(k)}$的检验。现在看一个案例，如何选择$k$来计算p值，如下所示：

1
2
3
4
5
6
7
8
alpha <- 0.05
i = seq(along=pvals)
mypar(1,2)
plot(i,sort(pvals))
abline(0,i/m*alpha)
##close-up
plot(i[1:15],sort(pvals)[1:15],main="Close-up")
abline(0,i/m*alpha)

1
2
3
k <- max( which( sort(pvals) < i/m*alpha) )
cutoff <- sort(pvals)[k]
cat("k =",k,"p-value cutoff=",cutoff)

结果如下所示：

1
2
> cat("k =",k,"p-value cutoff=",cutoff)
k = 11 p-value cutoff= 3.763357e-05

我们可以从数学上证明到这个程序可以将FDR控制在5%以下，具体的算法可以参考Benjamini-Hochberg在1995年的论文。这种新的程序计算的结果就是将原来得到的2个有统计学差异的数值提高到了11个。如果我们将FDR的值设为50%，那么这个数字会提高到1063。FWER无法提供这种灵活性，因为只要检测值的数量增加，都会造成FWER的值为1。

这里我们需要注意的是，在R中，已经有了计算FDR的函数，前面的那些复杂代码只是为了说明这种算法，在R中我们可以使用下面的代码进行计算，如下所示：

1
2
3
4
fdr <- p.adjust(pvals, method="fdr")
mypar(1,1)
plot(pvals,fdr,log="xy")
abline(h=alpha,v=cutoff) ##cutoff was computed above

我们也可以使用Monte Carlo模拟来确认一下FDR的值实际上小于0.05，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
alpha <- 0.05
B <- 1000 ##number of simulations. We should increase for more precision
res <- replicate(B,{
  controls <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
  treatments <- matrix(sample(population, N*m, replace=TRUE),nrow=m)
  treatments[which(!nullHypothesis),]<-treatments[which(!nullHypothesis),]+delta
  dat <- cbind(controls,treatments)
  pvals <- rowttests(dat,g)$p.value
  ##then the FDR
  calls <- p.adjust(pvals,method="fdr") < alpha
  R=sum(calls)
  Q=ifelse(R>0,sum(nullHypothesis & calls)/R,0)
  return(c(R,Q))
})
Qs <- res[2,]
mypar(1,1)
hist(Qs) ##Q is a random variable, this is its distribution

上述是Q值的直方图（假阳性除以显著性特征值的数目），现在看一下FDR值，如下所示：

1
2
FDR=mean(Qs)
print(FDR)

结果如下所示：

1
2
3
> FDR=mean(Qs)
> print(FDR)
[1] 0.03813818

这个FDR的值小于0.05，这个结果是符合预期的，因为我们需要保守一点，从而确保任何值的$m_{0}$的FDR都小于0.05，例如当每个假设检验都是零的极端情况，例如$m=m_{0}$。如果你重新模拟上述情况，你会发现FDR会增加。

我们需要注意下面的值：

1
2
3
> Rs <- res[1,]
> mean(Rs==0)*100
[1] 0.7

在模拟中，这个比例是0.7%，我们没有调用任何有显著性差异的基因。

在R中，p.adjust()函数可以选择一些算法来控制FDR，如下所示：

1
2
> p.adjust.methods
[1] "holm"       "hochberg"   "hommel"     "bonferroni" "BH"         "BY"         "fdr"        "none"

上面的这些方法不仅仅是估计错误率的不同方法，并且它们的估计值也不同，也就是说FWER与FDR不同。

直接FDR与q值

这里我们先回顾一下由John D.Storey在J. R. Statist. Soc. B(2002)中提到的结果。Storey与Benjamini-Hochberg方法的不同之处在于，前者不再提前设定$\alpha$水平。因为在一些高通量实验中，我们感兴趣的仅仅是找到一个基因列表用于验证这些基因，我们会事先考虑将那些p值小于0.01的基因都进行验证。我们人随后会考虑估计的错误率。通常使用这些方法，我们会确保我们的$R>0$。在上述定义的FDR里，当$R=V=0$时，我们指定$Q=0$，因此我们可以计算FDR如下所示：

在Storey提出的方法里，我们需要构建一个非空列表，也就是说$R>0$，那么我们计算阳性FDR(positive FDR)的公式如下所示：

第二点不同之处在于，Benjamini和Hochberg的程度是在最差的情况下进行控制，最差的情况是指零假设都为真（$m=m_{0}$的情况），Storey的方法则是让我们从数据中估计$m_{0}$。因为在高通量实验中，我们已经获得了如此多的数据，使Storey的算法成为了可能。这种算法的大致思路就是，挑选一个相对高值的p值cut-off，将它称为$\lambda$，并且假设那些p值大于$\Lambda$的检验在多数情况下其零假设是成立的，因此我们可以计算出估计值$\pi_0 = m_0/m$为：

还有其它比这种算法更复杂的算法，但是它们的思路都是一样的，例如当我们设定$\lambda=0.1$时，我们计算一下p值，如下所示：

1
2
3
4
5
hist(pvals,breaks=seq(0,1,0.05),freq=FALSE)
lambda = 0.1
pi0=sum(pvals> lambda) /((1-lambda)*m)
abline(h= pi0)
print(pi0) ##this is close to the trye pi0=0.9

1
2
> print(pi0) ##this is close to the trye pi0=0.9
[1] 0.9311111

根据这个估计，我们可以改变一下Benjamini-Hochberg程序，例如选择一个$k$作为最大值（$k$这里是最大的$i$），因此就有如下公式：

我们还有一种方法可以计算校正后的p值，即对每个检验计算q值(q-value)，如果一个特征计算的p值为$p$，那么q值就是一系列含有最小p值为$p$的特征值的估计pFDR。

除此之外，我们还有一种方法可以计算校正后的p值，即对每个检验计算q值(q-value)，如果一个特征最终计算的p值为$p$，那么q值就是一系列含有尽可能最小与$p$相等的p值的特征值的估计pFDR（原文是：However, instead of doing this, we compute a q-value for each test. If a feature resulted in a p-value of $p$, the q-value is the estimated pFDR for a list of all the features with a p-value at least as small as $p$.）。

上面这段我也不理解，后来查中文资料，根据Benjamini-Hochberg的算法，q值的定义如下所示：

对于m个独立的假设检验，它们对就的p值为：$p_{i},i=1,2,\cdots,m$

  1. 按照升序的方法对这些p值进行排序，得到：
  1. 对于给定的统计学检验水平$\alpha in(0,1]$，找到最大的$k$，使得：
  1. 对于排序靠前的$k$个假设检验，认为它们是真阳性，看下面的案例：

现在我们做了6个基因的检验，它们的p值如下所示：

Gene p-value
G1 p1=0.053
G2 p2=0.001
G3 p3=0.045
G4 p4=0.03
G5 p5=0.02
G6 p6=0.01

现在取检测水平$\alpha=0.05$，先把p值按从小到大的升序排列：

Gene p-value
G2 p2=0.001
G6 p6=0.01
G5 p5=0.02
G4 p4=0.03
G3 p3=0.045
G1 p1=0.053

取检测水平$\alpha=0.05$，现在我们找到一个值$k$，这个$k$满足以下公式：

当$k=1$时，$p_{(1)}=0.001 < 0.05*1/6=0.008333333$

当$k=2$时，$p_{(2)}=0.01<0.05*2/6=0.01666667$

当$k=3$时，$p_{(3)}=0.02<0.05*3/6=0.025$

当$k=4$时，$p_{(4)}=0.03<0.05*4/6=0.03333333$

当$k=5$时，$p_{(5)}=0.045>0.05*5/6=0.04166667$

当$k=6$时，$p_{(6)}=0.053>0.05*6/6=0.05$

从上面的计算过程可以知道，这个$k$等于4，也就是说，在FDR<0.05的情况下，G2，G6，G5，G4有差异。

到这里，只是说明了，G2，G6，G5，G4是有差异的，但是q值还没有算出来，继续计算：

前面我们已经把原始的p值按照从小到大的顺序排列好了，也就是[1] 0.001 0.010 0.020 0.030 0.045 0.053，其中最大的p值是0.053，它校正后还是这个值，也就是说这个值是校正后的最大p值，次大的p值是0.045，这个值需要校正，它排序是第5，那么校正的公式就是所有p值的数目（一共是6个p值）除以秩（这里是5），再乘以p值大小，也就是0.045*6/5=0.054，但是，这个值已经大于原来最大的p值（0.053）了，因此这个把它校正后的值也作为0.053，再看倒数第3个值，即0.03的校正p值，即0.03*6/4=0.045，它小于0.053，因此它校正后可以是0.045，现在依次校正剩下的值：

0.02*6/3=0.04，0.01*6/2=0.03，0.001*6/1=0.006，也就是说校正后的p值（就是q值），按从小到大的顺序排列就是：0.006，0.03，0.04，0.045，0.053，0.053，从结果中我们可以发现，前4个是有差异的，它们都小于0.05，也就是说G2，G6，G5，G4有差异。

其实公式就是：

以上是BH法进行q值计算的过程，R中可以使用p.adjust()函数计算q值，所示：

1
2
3
4
5
6
7
8
9
10
11
12
p1 <- 0.053
p2 <- 0.001
p3 <- 0.045
p4 <- 0.03
p5 <- 0.02
p6 <- 0.01
p0 <- c(p1,p2,p3,p4,p5,p6)
p.adjust(p0,method = "BH")
p.adjust(sort(p0),method = "BH")
p.adjust(sort(p0),method = "fdr")
sort(p.adjust(p0,method = "BH"))

结果如下所示：

1
2
3
4
5
6
7
8
> p.adjust(p0,method = "BH")
[1] 0.053 0.006 0.053 0.045 0.040 0.030
> p.adjust(sort(p0),method = "BH")
[1] 0.006 0.030 0.040 0.045 0.053 0.053
> p.adjust(sort(p0),method = "fdr")
[1] 0.006 0.030 0.040 0.045 0.053 0.053
> sort(p.adjust(p0,method = "BH"))
[1] 0.006 0.030 0.040 0.045 0.053 0.053

从上面的计算结果看来：

  1. method="BH"等于method="fdr"，结果没有区别；
  2. 在使用p.adjust()函数计算q值时，不用对原始数据进行排序，如果想要结果按从小到大的序列排列，那么就排序原始值，或计算后的值均可。

回到原文。

在R中，qvalue包中的qvalue()函数可以用来计算q值，如下所示：

1
2
3
4
library(qvalue)
res <- qvalue(pvals)
qvals <- res$qvalues
plot(pvals,qvals)

估计一下$\hat{\pi_{0}}$，如下所示：

1
2
> res$pi0
[1] 0.8813727

练习

P274

基础探索数据分析

与低通量数据相比，高通量数据的一个被低估的优点就是它很容易展现数据。例如当我们有了海量的高通量数据后，我们很容易发现那些在少量数据时并不明显的问题。研究这些数据的一个强有力工具就是探索性数据分析(EDA,exploratory data analysis)。这一部分我们就来了解这方面的内容，可以参考作者的Github上的原文。

火山图

我们使用前面数据做的t检验结果来看一下火山图，如下所示：

1
2
3
4
5
6
library(genefilter)
library(GSE5859Subset)
data(GSE5859Subset)
g <- factor(sampleInfo$group)
results <- rowttests(geneExpression,g)
pvals <- results$p.value

我们还可以从一个数据集中生成一些p值，这些数据集中的零假设为真，如下所示：

1
2
3
4
m <- nrow(geneExpression)
n <- ncol(geneExpression)
randomData <- matrix(rnorm(n*m),m,n)
nullpvals <- rowttests(randomData,g)$p.value

我们前面提到过，当我们报告效应大小(effect size)时，如果仅仅计算p值则会造成一些错误。我们可以通过画一个火山图来展示高通量数据的统计结果。火山图背后的思想是，它能展示出所有的特征值（这里指的是基因表达谱）。火山图的y轴是-log(base 10) p-value，x轴是效应大小(effect size)。当我们经过-log(base 10)转换后，那些有着极显著差异的特征值就会出现在火山图的上方。这里使用log转换是因为，我们可用log转换可以更好地区分一些非常小的数据，例如区分0.01和$10^6$，此时我们来绘制一个火山图，如下所示：

1
2
plot(results$dm,-log10(results$p.value),
xlab="Effect size",ylab="- log (base 10) p-values")

p值直方图

另外一种看整体数据的思路就是绘制p值的直方图。当我们的数据完全无效时，那么p值的直方图是服从均匀分布的，而我们的数据有效时，我们可以发现较小的p值频率较高，如下所示：

1
2
3
4
library(rafalib)
mypar(1,2)
hist(nullpvals,ylim=c(0,1400))
hist(pvals,ylim=c(0,1400))

左侧的p值是无效数据产生的p值，它服从均匀分布，右侧的p值则是则是我们基因表达谱的数据。

当我们预期大多数假设都是无效时，我们就不会看到服从均匀分布的p值，它也许是一些异常属性的指标，例如相关的样本。如果我们对结果时行置换检验，并计算出p值后，如果样本是独立的，那么我们应该会看到均匀分布，但是，我们的数据却看不到这个结果：

1
2
3
permg <- sample(g)
permresults <- rowttests(geneExpression,permg)
hist(permresults$p.value)

在后面部分中，我们会了解到这个数据集中的每个检验并不是独立的，因此这里用于检验的假设（我们使用的是t检验，而t检验的前提就是样本独立）是不正确的。

箱线图与直方图

高通量数据实验中，我们会检测每个实验单元的数千个特征值，我们前面也提到了，使用箱线图可以辅助我们来判断这些数据的质量。例如，如果一个样本的分布完全不同于剩余的样本，那么我们就可以认为，这个样本存在着一定问题。虽然一个完全的完全的变化可能是由于真正的生物学差异引起的，但是多数情况下，这就是技术因素造成的。这里我们使用Bioconductor中的基因表达数据，为了模拟出一组异常的数据，我们会对其中的一个样本进行log2转换，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
#BiocManager::install("Biobase")
#devtools::install_github("genomicsclass/GSE5859")
library(Biobase)
library(genefilter)
load("GSE5859.rda")
data(GSE5859)
ge <- exprs(e) ##ge for gene expression
ge[,49] <- ge[,49]/log2(exp(1)) ##immitate error
library(rafalib)
mypar(1,1)
boxplot(ge,range=0,names=1:ncol(e),col=ifelse(1:ncol(ge)==49,1,2))

运行过重中发现，GSE5859这个包无法安装，因此可以去官网下载原始文件，直接加载到RStudio中，另外在加载data(GSE5859)时会出错，不用管，运行就行，图形如下所示：

从上图我们可以看到，样本数据大，很难看清楚中间的箱子形状，但是我们很容易发现有一个样本异常，除此之外，我们可以用另外一种方式来展示一下数据，这个数据不画箱子，如下所示：

1
2
qs <- t(apply(ge,2,quantile,prob=c(0.05,0.25,0.5,0.75,0.95)))
matplot(qs,type="l",lty=1)

这种图形可以称为kaboxplot，因为这是由Karl Broman首次使用的，它绘制的是0.05，0.25，0.5，0.75和0.95分位数的图形。

我们也可以绘制直方图。因为我们的数据很多，因此可以使用很窄的间隔(bin)与柱子，然后将这些柱子进行平滑处理，最终绘制成平滑直方图(smooth histogram)。我们重新再校正这些平滑曲线的高度，那么在$x_{0}$处的高度与基本单元构成的面积内，它们的点的数目就都比较接近，如下所示：

元区域内的点的数目就与这个基本单元的面积接近积，如下所示：

1
2
mypar(1,1)
shist(ge,unit=0.5)

MA图

散点图(scatterplot)与相关性(correlation)不是检测重复性问题的最佳工具。检测重复性更好的工作就是检测两次之间的差值，如果重复性好，那么这些差值应该是一样的。因此，一种更好的绘图工具就是将散点图旋转一下，这个散点图的y轴上差值，x轴是平均值。这种图最初被叫做Bland-Altman图，但在遗传学中它经常被称为MA图。MA的名称来源于图形的内容，M表示减(minus)，表示两值之差（有的时候是log值之差），A表示平均(Average)，现在绘制一下MA图

1
2
3
4
5
x <- ge[,1]
y <- ge[,2]
mypar(1,2)
plot(x,y)
plot((x+y)/2,x-y)

左图是常规的相关图，右图是MA图，需要注意的是，当我们把左图进行旋转后，这些数据的差异的SD就非常直观了：

1
2
sd(y-x)
[1] 0.2025465

左图的散点图显示这两个样本的相关性很强，但是显示的信息有限。

参考资料

  1. 如何理解与计算FDR？（第二版）